RKKY Interaction in Graphene from Lattice Green's Function 
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We study the exchange interaction J between two magnetic impurities in graphene (the RKKY 
interaction) by directly computing the lattice Green's function for the tight-binding band structure 
for the honeycomb lattice. The method allows us to compute J numerically for much larger distances 
than can be handled by finite-lattice calculations as well as for small distances. In addition, we 
rederive the analytical long-distance behavior of J for linearly dispersive bands and find corrections 
to the oscillatory factor that were previously missed in the literature. The main features of the 
RKKY interaction in graphene are that unlike the J oc (2IcfR)~ 2 sin(2fcF-R) behavior of an ordinary 
2D metal in the long-distance limit, J in graphene falls off as 1/R 3 , shows the 1 + cos ((K — 
K').R)-type oscillations with additional phase factors depending on the direction, and exhibits a 
ferromagnetic interaction for moments on the same sublattice and an antiferromagnetic interaction 
for moments on the opposite sublattices as required by particle-hole symmetry. The computed J 
with the full band structure agrees with our analytical results in the long-distance limit including 
the oscillatory factors with the additional phases. 

PACS numbers: 75.30.Hx; 75.10.Lp; 75.20.Hr 



I. INTRODUCTION 

Graphene has attracted considerable attention recently 
due to its linear energy dispersion, where the excitations 
are massless Dirac Fermions, which could lead to phys- 
ical behavior different from that of the standard two- 
dimensional systems. The Ruderman-Kittel-Kasuya- 
Yosida (RKKY) interaction is the exchange interaction 
between two magnetic impurities mediated by the con- 
duction electrons of the host and is a fundamental quan- 
tity of interest J^- Earlier works on the RKKY interaction 
in graphene 3 ""™ have mostly used a continuum model with 
a linearly dispersive band structure with two Dirac cones 
at the corners of the Brillouin zone or have used exact 
diagonalization on a finite size latticed As pointed out 
by Saremi)2>£ for a bipartite lattice with nearest-neighbor 
interactions, particle-hole symmetry leads to a ferromag- 
netic interaction on the same sublattice and an antifer- 
romagnetic interaction for the opposite sublattices. 

The continuum limit, linearly-dispersive Dirac cone ap- 
proximation is expected to be good in the long-distance 
limit R — ¥ oo and allows an analytical solution. However, 
it requires the use of a cutoff function, without which 
the contributions from the higher energy states add up 
to produce a diverging and oscillatory result. If a sharp 
energy cutoff is used, it produces a J that oscillates with 
distance violating the particle-hole symmetry. To address 
this problem, a cutoff function approach 3 - has been used, 
where the higher-momentum contributions are damped 
out slowly, with the length scale of damping taken to in- 
finity as the limiting case. While this approach yields a 
reasonable result, it is not a priori obvious if some sys- 
tematic error is not introduced by such a procedure. In 
fact, using exact diagonalization on finite lattices, Black- 
Schaffer 7 extracted J values which differed from Saremi's 
results both in the prefactors and, for Jab , in the oscilla- 
tory factor as well. On the other hand, the finite lattice 



calculations in turn suffer from the deficiency that the 
distance between the impurities can not be too large and 
extra interactions between the moments get introduced 
due to the supercell geometry. 

In order to address these issues, we use a direct com- 
putation of the lattice Green's functions to compute the 
RKKY interaction with the full tight-binding band struc- 
ture. The method allows us to numerically calculate the 
RKKY interaction for very large distances, which is im- 
practical to obtain from the finite lattice calculations. 
In addition, we obtain analytical results for the long- 
distance behavior of J by using an approach slightly dif- 
ferent from that of Saremi, 3 which allows us to obtain the 
proper phase factors in the RKKY oscillations that were 
missed in the previous works. We note that for obtain- 
ing the proper oscillatory factors of </, it is important to 
include carefully the phase factors of the electronic wave 
functions around the Dirac points. 



II. FORMULATION 

We consider the tight-binding Hamiltonian for 
graphene including a contact interaction with the mag- 
netic centers, viz., 



H 



tij C ia C j<T 



h.c. — A S p 



(1) 



where i is the combined site-sublattice index, (ij) de- 
notes summation over distinct pairs of nearest-neighbor 
sites, a denotes the electron spin, the p summation runs 
over the magnetic centers, and s p = ( fr/2) Ylnv c \^^Cpv 
is the itinerant electron spin density. The results can 
be easily generalized if the hopping integrals tij are re- 
tained beyond the nearest neighbor. However, direct nu- 
merical computations showed that neglecting the higher 
neighbor terms does not change J significantly, since the 



FIG. 1: The Hamiltonian near the Dirac 

points in the first Brillouin zone with <f>(q) = 

_ e i(7r/3-9,) ) gi(7r/3+fl,) g-ifl, _ e i(2n / 3+6 q ) ^ gi(2ir/3-9») 

for the first through the sixth Dirac points (ffi to i^) as 
labelled in the figure. The small momentum q is the deviation 
from the corresponding Dirac point (k — Kd + q) and 6 q is 
the polar angle of q with respect to Ki — Ki chosen as the 
x direction as shown in the figure. These phase factors are 
important as they determine the oscillatory behavior of Jab 
as discussed in the text. 

magnitudes of the hopping beyond the nearest-neighbors 
are relatively small in graphene £ The nearest-neighbor 
hopping parameter in graphene is t — —2.56 eV as ob- 
tained by fitting the tight-binding bands to the density- 
functional band calculations.— 

In the basis of the Bloch functions {c\ ao = 
N~ 1/2 J2i e lfe - (fll+Ta) cJ Q(T ) of the two sublattices, a = A 
or B, the unperturbed Hamiltonian is given by 

Hk = ( /*(*) / (f ) ) ' (2) 

where f(k) = t (e ifc ' dl + e lk ' d2 + e ifc ' d3 ). The Hamil- 
tonian expression near a Dirac point takes the form 
f(k) = vpq <f>(q), where 4>(q) is a phase, different at dif- 
ferent Dirac points as indicated in Fig. [TJ and q is the 
deviation of the Bloch momentum k from the neighboring 
Dirac point, q = k — Kf>. Diagonalizing the Hamiltonian, 
one finds a linear dispersion near the Dirac points, viz., 
E q = ±\vf\q with the Fermi velocity vf = Sat/2, which 
is defined to be negative throughout this paper, t being 
negative, and here a is the bond length. 

In the linear response theory, the exchange interac- 
tion may be obtained by first computing the perturbed 
wave functions due to the magnetic impurity located 
at the origin from the Lippmann-Schwingcr equation 
I*) = |*°) + GV\^), from which the energy due to 
the second magnetic impurity located at R is com- 
puted from the first-order perturbation theory, so that 
E(R) = {^\H int (R)\9}. For a contact interaction be- 
tween the magnetic impurity and the conduction elec- 
trons, Hint = —A S ■ s5(f), E(R) may be written in the 



FIG. 2: The graphene honeycomb lattice with two different 
sublattices, shown as full and open circles. The figure shows 
the orientation of the Brillouin zone and two common direc- 
tions in the direct lattice (zigzag and armchair). 

Heisenberg form 

E{R) = JSx ■ § 2 , (3) 

where the exchange integral J = (\ 2 h 2 /4) %(0,i?). The 
susceptibility x( r ? r ') — 5n(r)/SV(r') is written in terms 
of the unperturbed retarded Green's functions 

X (r,r') = ~ J F dEIm[G°(r,r\E)G°{r l ,r,E)}. (4) 

Here G° is the Green's function for a single spin chan- 
nel and is spin independent, while in the definition of 
the susceptibility, 6V(r') is a spin-independent perturba- 
tion and 5n(r) is the change in the total charge density 
including both spins. 

This result can be easily extended to the case of 
graphene to yield 

Xa/s(r, r') = -1 p dE Im[G^(r, /, E)G$ a (r' , r, E)], 

where a,/3 are the sublattice indices, A or B, r, r' de- 
note the lattice positions of the two magnetic centers 
located on the sublattice a and /3, respectively, and 
the sublattice susceptibility is as usual Xapir^r') = 
Sn a (r) I '6Vp(r'). The expression for the susceptibility 
is obtained by noting that the charge density is given 
by n a (r) = — — J Ef dE Im G QQ (r, r, E) and obtaining 
the charge difference Sn a (r) induced by the perturbation 
SV p (r') from the Dyson equation G = G° + G°VG. The 
exchange interaction J a p(R) between impurities located 
at the sites (a, 0) and ((3, R) is then given by 

J al }(R) = ^-x a p(0,R)- (6) 

The calculation of the exchange interaction thus boils 
down to the computation of the lattice Green's functions 
and a quadrature over the energy following Eq. [5J 
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III. CALCULATION OF THE GREEN'S 
FUNCTIONS 

We compute the real-space Green's function by two 
different approaches, viz., the direct integration method 
and the recursive technique of Horiguchiii£ In the first 
method, the real-space Green's functions are calculated 
by numerically integrating the momentum-space Green's 
functions 



G%(r,r',E) 



1 







BZ 



where 



G°{k,E) 



^*BA ^BB 

E + irj + n k 
(E + ir,)*-f*(k)f(ky 



(E + i V - H k 



(8) 



The Brillouin zone integral in Eq. [7] was evaluated 
by taking up to = 10 6 fc-points in the full Brillouin 
zone and a small value for the infinitesimal parameter 
7] = 0.005 is used. This method is straightforward and 
computationally robust but slow, while the Horiguchi re- 
cursive technique is fast, but it has stability problems 11 
for larger distances. 

It is worth noting the symmetry properties of the 
Green's functions, which immediately follow from the 
above equations and the expression for "H/c, viz., that 
G AA (R,0,E) = G BB (R,0,E) and G AB (0,R,E) = 
Gba(R,0,E), leading to the results, which is also ob- 
vious on physical grounds: 

Jaa(R) = Jbb(R) and J AB (R) = Jba(-R)- (9) 

In the second method, the Horiguchi recursive 
technique^ the Green's functions for the honeycomb 
lattice is expressed in terms of those for the triangular 
lattice, which in turn are expressed in terms of the ellip- 
tic integrals. For example, the expression for the on-site 
Green's function is given by 



lA (0,0,z) = — g(z')K(k(z% 



(10) 



where z = E+irj, z' = (z 2 -3)/2, g(z') = 8((2z' + 3) 1 / 2 - 
l)- 3 / 2 ((2z' + 3) 1 / 2 + 3)- 1 / 2 , k(z') = 2~ 1 g(z')( 2z' + 3 ) 1 / 4 , 
K(k) = K(k) for Im k < and K(k) + 2iK(sJ\ - k 2 ) for 
Im k > 0, and K(k) is the elliptic integral 



K{k) 



7r/2 



(1-fc 2 sin 2 0) 1 / 2 



(10. 



(11) 



The elliptic integral with complex modulus was evaluated 
following established procedures and the Arithmetic- 
Geometric Mean Metho d 12 ' 13 and r\ ~ 10~ 6 was used. 

The computed result for the on-site Green's function 
obtained from Eq. [TU] is shown in Fig. [3J and one 
sees the familiar density-of-states, which is proportional 
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FIG. 3: On-site Green's function Re G° AA {0, 0, E) (dashed 
line) and — Im G AA (0, 0, E) (full line) computed from Eq. [TUJ 
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FIG. 4: (Color online) Real (dashed line) and Imaginary (full 
line) parts of the Green's function G AA (R,0,E) with R = 
7\/3a(l,0) obtained from the Horiguchi recursive method. 
Red line shows the product of the Real and Imaginary parts, 
which is displaced by +0.25 along the y-axis. 



to the imaginary part. Similarly, one can compute the 
Green's functions for a few lattice vectors R (specifically, 
R = (l,m) = (0,0), (2,0), and (4,0), for the triangular 
lattice), from which the remaining Green's functions for 
both the triangular as well as the honeycomb lattice can 
be computed using the recursion relations. 

In Fig. HI we have plotted the Green's function 
G° AA {R,0,E) with a specific R = 7^(1,0) along the 
zigzag direction and also the product of the real and the 
imaginary parts, which enters as the integrand in the cal- 
culation of J (Eq. [5]) , to be integrated over the occupied 
states between E = —3 and zero. As seen from the fig- 
ure, the integrand is a rapidly oscillating function, with 
a small net result for J. 



IV. EXCHANGE INTERACTION 

Before we present the results of the full calculations, 
we derive the long-distance behavior of J with the lin- 
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early dispersive band structure. This is reasonable as the 
long-distance behavior is necessarily determined by the 
small momentum states, for which the linear-band ap- 
proximation is excellent. As has been pointed out in the 
literature*^ the RKKY interaction shows oscillations as 
a function of distance because of the interference between 
charge densities originating from the two Dirac points in 
the Brillouin zone. However, there is no consensus re- 
garding the form of the oscillation and we find important 
differences from the earlier results. 



A. Magnetic impurities on the same sublattice 

To derive the long-distance oscillatory behavior of the 
RKKY interaction, we first obtain the Green's functions 
for small energies (which necessarily determine the long 
distance behavior) and use these to evaluate the integral 
in the expression Eq. [3] for the susceptibility. For small 
energies, the contribution to the integral in Eq. [7] comes 
from the two Dirac points in the BZ, so that the equation 
becomes 



C'sffi, 0, E) 



1 



9. 



BZ 



d 2 q e^ R 



x W K - R G%(q + K,E) + e iK ' R G° ap ( q + K' , 25)£L2) 

where K and K' are the two Dirac points and R is the 
distance between the two magnetic centers. 

For a — (3, i.e., if the two sublattices are the same, 
these two Green's functions in the kernel become the 
same, viz., G° AA (q,E) = (E + irj){{E + iij) 2 - vpq 2 )- 1 
as seen from Eq. [8] The exponential factor may be 
expanded (Jacobi-Anger expansion^) in terms of the 
Bessel's functions 



Jq-R 



J (qR) + 2 V i n J„(qR) cos[n(9 q - 6 a )] (13) 



n=l 



and then integrating over 8 qi we get the result 

G° AA (R, 0, E) = (e lK - R + e lK ' R )g AA (R, E), (14) 
where 

2tt 



9aa(R, E) 



BZ Jo 



qdqJ (qR)G AA (q,E) (15) 



and a momentum cutoff q c has been introduced. We get 
a similar expression for G° AA (0, R, E) and plugging these 
in Eq. [5J we immediately get 

X aa(0, R) = Iaa(R) x (1 + cos[(K - K') ■ R\) (16) 
with the prefactor 



Iaa(R) 



dE Im [g AA (R,E)f 



(17) 



The large distance behavior of J is controlled by 
small momentum states and one may try to evaluate 



this by taking the cutoff q c — >• oo for the ease of per- 
forming the integrals. In this case, the integral in Eq. 
[TBI can be expressed in terms of the modified Bessel's 
function of the second kind, viz., g AA — —2ttQ, 



BZ 



E/vp x K {iER/vF)^^ This in turn may be ex- 
pressed in terms of Bessel and Neumann functions 17 with 
real arguments to yield the kernel Im [g AA (R, E)] 2 — 
(2TT 4 )n B 2 z v F 2 R- 2 y 2 J Q (y)Y {y), where y = ERv F \ Thus 
Eq. [T"7l becomes 



Iaa{R) = 



8n J 



-R- 



dyy 2 J (y)Y (y). (18) 



r BZ v F 

The R~ 3 dependence clearly emerges; however, the 
integral does not converge. Following Saremi^ we 
multiply the integrand by a cutoff function f(y,yo), 
perform the integral, and take the limit j/o ~^ 
oo. We have tried three different cutoff functions, 
viz., f(y,ya) = exp (-y/y ), exp (-y 2 /iJo), or 
VoiVo + y 2 )" 1 : an d in each case find the same limit: 
lim^o-KM J Q °° dyy 2 J (y)Y (y)f(y,y ) = 1/16. Thus we 
immediately get Iaa(R) — 9(647r£)~ 1 a 3 /i? 3 , which leads 
to the exchange interaction for the same sublattice 



Jaa = ~C X 



1 + cos((K - K') ■ R) 

(RH 3 : 



(19) 



where C = — 9A 2 ft 2 /(2567rf) is a positive quantity, since 
t < in graphene. We note that the oscillatory factor 
is identical to the expression derived by Saremij^ how- 
ever, we have an additional scaling factor of 3/2 for the 
magnitude of Jaa- From finite-size calculations, Black- 
Schaffer— extracted a scaling factor different from ours 
or that of Saremi; however, we note that such factors 
are difficult to extract from numerical results, especially 
from finite-size calculations. 

The expression for J A a, Eq. [TJU is valid for all di- 
rections including the zigzag and the armchair directions 
and K and K' are any two adjacent Dirac points in the 
Brillouin zone. It is easy to see that while the oscilla- 
tory factor 1 + cos((K — K') ■ R) repeats in triplets as 2, 
1/2, 1/2, ... with distance R along the zigzag direction, 
it is always two for the armchair direction. Because of 
this, the magnitude of Jaa oscillates for the zigzag direc- 
tion but not for the armchair direction, always however 
remaining ferromagnetic as required by the particle-hole 
symmetry. 

The calculated results for Jaa using the full band 
structure following methods of Sees. II and III are shown 
in Fig. [5] for the zigzag direction and in Fig. [6] for the 
armchair direction. Both follow the long-distance behav- 
ior of Eq. [19] quite well beginning with surprisingly small 
distances. 



B. Magnetic impurities on two different sublattices 

We now turn to Jab , where the two impurities are lo- 
cated on different sublattices. The needed Green's func- 
tions, in the small q limit, are obtained from Eq. [5] to 
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100 150 200 250 
R/a 



FIG. 5: (Color online) Exchange interaction Jaa between two 
impurities on the same sublattice. Black lines are the results 
with the full tight-binding band structure, while the red lines 
indicate the long-distance behavior as obtained from Eq. [T5] 
using linear dispersion. The inset shows the log plot showing 
the long-distance R~ 3 behavior, while there are noticeable 
differences for small R, especially visible in the inset. Note 
that since t is negative for graphene, Jaa is also negative 
indicating a ferromagnetic interaction. 
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FIG. 6: (Color online) Same as Fig. [5]for the armchair direc- 
tion. For this direction, consistent with Eq. I19I there are no 
oscillations in J unlike the zigzag direction (Fig. [5j. 



yield 



G° BA (q + K D ,E) = ± 



-»(7r/3±e„) 



(E + ir]) 2 - v 2 F q 2 



(20) 



where ± signs are for the two Dirac points Kp = K and 
K' 7 respectively. Here, we have chosen a Brillouin zone 
that includes the K \ and K2 points and the correspond- 
ing phase factors in the Hamiltonian have been retained 
(see Fig. [T]). The next step is to obtain G BA (R,0,E) 
by the momentum space integration using Eq. I12I The 
same Jacobi- Anger expansion for e lq ' R (Eq. fT5|) may be 
used as before except that now the extra phase factor 
e ±i8 q a pp ears i n the angle integral while performing the 
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FIG. 7: (Color online) Jba for the zigzag direction. Black 
solid line is the result of the full calculation, while the red 
dashed line is obtained from Eq. I26I 
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FIG. 8: (Color online) Same as Fig. [7]for the armchair direc- 
tion. 

integration in Eq. Q21 Using the result 

'0 if n 7^1, 




d0 q e 



" cos[n(6 q - Or)} = 



ne ±WR if n = 1 



(21) 



where Or is the polar angle of the distance vector 
R as defined in Fig. [21 we get after some algebra, 
the result: G BA (R,0,E) — a gBA{R,E) 7 where a = 



r/3f e i(K-R-8 R ) _ e i(K'-R+e R )\ an( j 

g B A{R,E) = — / dq 



n 



BZ Jo 



(E + i-q) 2 — v 2 F q 2 ' 



(22) 



Similarly we find G° AB (0,R,E) = ~a* g B A(R,E). Fi- 
nally, using Eq. O the susceptibility becomes 

Xba{R, 0) - Iba(R) x (1 + cos[(K — K') ■ R + tt — 26 R \) 

(23) 

with the prefactor 



4 f E " 

Iba{R) = -J dElm [g BA {R, E)Y 



(24) 



We can now proceed to evaluate Jba in the long- 
distance limit in a fashion similar to the previous subsec- 
tion. We find that qba = 27rf2^ z x E/v 2 F x Ki(iER/vf), 
where K\ is the first-order modified Bessel function of the 
second kind. Expressing this in terms of Bessel and Neu- 
mann functions and using a cutoff function as before, we 
finally get 



Iba(R) 



n 2 BZ v F 



! x lim / dyy 2 Ji(y)Yi(y)f(y,y ) 

(25) 

The result for this integral is -3/16, so that collecting all 
terms, the exchange interaction becomes 



Jba = 3C x 



1 + cos((K - K') ■ R + 7T - 29 R ) 



(26) 



Note that the equation is valid for any direction of R 
and for any choice of the two Dirac points K and K 1 
(they have to be adjacent to each other of course, so 
that they are within a single unit cell in the reciprocal 
lattice), so long as we define the angle Or to be with 
respect to the chosen K 1 — K vector. One can check that 
these equations yield the same results for J for different 
equivalent directions R as expected from symmetry. 

Note that Jba has the extra factor tt — 29r in the 
argument of the cosine as compared to Jaa, which comes 
from the interference of the contributions to the Green's 
functions Eq. Q2] from the two Dirac points. For the 
zigzag direction along K' — K, the oscillatory part of 
Eq. [25] agrees with the B lack- S chaffer result^ since the 
angle Or vanishes for large R. However, our result is 
valid for all directions and furthermore for the armchair 
directions, the angle 9r is never zero (see Fig. [2]), so that 
this phase factor must be retained in Eq. [211 

The oscillatory behaviors seen in the full calculations 
for Jab presented in Figs. [7] and are contained in 
Eq. [26] For the zigzag direction, taking K — K' — 
47r(3v / 3a)~ 1 (— 1, 0) and R = d% + V3anx, where n is 
an integer, the oscillatory factor in Eq. [23] becomes 
1 + cos[(4n - l)7r/3 + 26 R ]. In the limit R ->• oo, 9 R 0, 
so that this factor repeats in the sequence of the triplet 
numbers: 0, 3/2, 3/2 as R is increased. For a finite R, 
Or 7^ 0, so that we never get exactly the zero in the 
triplet, but rather a small number, which is faithfully re- 
produced in the full calculations shown in Fig. [7] where 
two values of J are close in magnitude, while the next 
one is lower by about three orders of magnitude. For the 
armchair direction, R = 2~ 1 a(3n + l)(\/3, 1), n being 
an integer and 9r — ir/6, so that the oscillatory factor 
in Eq. [23] is always two. The exchange interaction thus 
changes smoothly with distance without any oscillations 
as seen from Fig. [8] 



C. Interaction between impurities on plaquettes 

The plaquette impurities, where the magnetic impu- 
rities are located at high-symmetry points rather than 
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FIG. 9: Magnetic interaction J for the hexagonal plaquette 
along the armchair direction. 
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FIG. 10: Magnetic interaction J for the hexagonal plaquette 
along the zigzag direction. 



at single lattice sites, are of interest because a number 
of atoms and small molecules may be favored to occupy 
such positions. For instance, the impurity may be located 
midway at a bond center and may interact with the two 
neighboring sites on the bond. In these situations, the 
interaction term in Eq. [TJ may be replaced by 



Hi 



(27) 



where the summations are performed over the lattice sites 
with which the impurity spins Sj interact (six sites if the 
impurity is located at the hexagon center and two, if it 
is on the bond center). Assuming that the interaction 
strength is small (A << |t|), the interaction J for these 
plaquette impurities may be obtained by simply summing 
over the individual site-site interactions already obtained 
in Sec. Ill, so that J = J2 P , P > Jpp'- 

These results are shown in Figs. [H]and[TU]for the hexag- 
onal plaquettes and in Fig. [TT] for impurities located on 
the bond centers. For the hexagonal plaquettes, as has 
been pointed out earlier^ the oscillating cosine factors 
present in the site interactions, J aa and Jab, cancel out 
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FIG. 11: (Color online) Magnetic interaction Jbond between 
two impurities located on bond centers obtained from the full 
calculation (solid line) and the linear-band long-distance limit 
Eq. [29] (red dots). 

both for the armchair and the zigzag hexagonal plaque- 
tte, leading to the single result valid for both cases: 

J P ia q = 36C x (R/a)- 3 , (28) 

which is a net antiferromagnetic interaction. 

For impurities located midway between the bond cen- 
ters along the zigzag direction, we find an oscillating in- 
teraction, in the long-distance limit, 

■Wfi)=4 C x '- 2c °^- g ''- R > . (29) 

(R/aY 

As distance R is increased, the numerator changes with 
the repeat sequence of -I, 2, and 2, leading to an in- 
teraction that is antiferromagnetic every third site and 
ferromagnetic otherwise. The results from the full calcu- 
lation is compared to the long-distance limit result, Eq. 
[29J in Fig. [TO 

V. SUMMARY AND DISCUSSIONS 

In summary, we studied the RKKY interaction be- 
tween the magnetic impurities in the honeycomb lat- 
tice by evaluating the Green's functions for the tight- 
binding Hamiltonian by the direct summation method, 
which worked well for all distances but is computation- 
ally slow, or the Horiguchi recursive technique, which is 
a fast method, but has stability problem for large R. For 
distances, where both methods worked, the results agreed 



with each other quite well. These methods are comple- 
mentary to the finite-lattice calculations^, however, the 
direct summation method allows the calculation of J for 
much larger distances with modest computational efforts. 
By carefully considering the phase factors of the wave 
functions around the Dirac cones, we have also obtained 
the analytical long-distance limits of J, Eqs. Q1J] and 121)1 
which are, although similar in form to previous results^i 
have important corrections in terms of additional phase 
factors in the oscillating term. All such phase factors 
were faithfully reproduced in our numerical calculations 
of J using the full tight-binding band structure. We 
found that the long-distance limit is reached for quite 
small distances, of the order of a few lattice constants. 

In addition to the nearest-neighbor model, we have also 
studied the effect of the further-neighbor electron hop- 
ping, but these produced negligible differences as might 
be expected, since the strengths of the higher-neighbor 
hoppings are quite small in graphene £ For the hexago- 
nal plaquette impurities, J is always antiferromagnetic, 
while for the bond impurities, the sign oscillates. Given 
that the magnitude of J falls off quite rapidly with dis- 
tance, nearest-neighbor J would dominate, so that spin 
chains based on hexagonal plaquette sites are prediced 
to be antiferromagnetic, while those based on the bond 
sites would be ferromagnetic. 

In trying to design an experimental system to observe 
the RKKY interaction, one must carefully select a proper 
system. At first sight, it might appear that a mag- 
netic adatom such as Co or Fe deposited on top of the 
graphene sheet would interact via the RKKY interaction. 
However, in addition to introducing the needed localized 
magnetic moments from the d electrons, the outermost 
s electrons of the adatom are transferred to the host, 
where they are added to either the conduction band or 
they may form weakly localized states around the adatom 
site. These extra electrons will modify the RKKY inter- 
action. The challenge is therefore to come up with a 
system, perhaps a simple molecule, that has a magnetic 
moment which interacts with the graphene lattice, but 
one that does not alter the electronic structure by con- 
tributing extra electrons to the graphene sheet. On the 
other hand, scaling arguments^ as well as renormaliza- 
tion group calculations^ indicate the lack of a Kondo 
effect below the critical coupling J c ps 3.5 eV, which is 
quite strong, so that the RKKY interaction should dom- 
inate. 
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